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Abstract. I review recent progress from iV-body simulations in our under- 
standing of the secular evolution of isolated disk galaxies. I describe some of the 
recent controversies in the field which have been commonly attributed to numer- 
ics. The numerical methods used are widely used in computational astronomy 
and the problems encountered are therefore of wider interest. 



1. Introduction 

The fragility of disk galaxies suggests that a significant fraction of their history 
was spent experiencing mild external perturbations. Thus their evolution is 
likely to have been driven partly by internal processes. Dissipational gas physics 
probably plays a prominent role in this evolution but the effects of collisionless 
evolution very near equilibrium are also significant. Collisionless secular evolu- 
tion is driven largely by non-axisymmetric structures, especially by bars, which 
represent large departures from axisymmetry, and are thus efficient agents of 
mass, energy and angular momentum r e distribution. Moreov er bars occur in 
over 50% of disk galaxies (lKnapenlll999l: lEskridge et alJl2000T ) and simulations 
starting from the 1970s dMiller et al.1 Il97fl: iHohl 197ll) have sh own that they 
form readily via global instabilities (|Kalnaislfl97l lToomrel lT981 ) . Here I review 
recent progress in understanding the secular evolution of isolated disk galaxies. 

2. Numerical Methods 

At present there is considerable debate on the role numerics play in the evolution 
of iV-body simulations. This review therefore treats collisionless isolated galaxy 
simulations as the necessary prerequisites for confidence in any other iV-body 
study of galaxy formation and evolution. The physics of isolated collisionless sys- 
tems consist only of Newton's law of gravity and Newtonian mechanics, both of 
which are well understood. Other than initial conditions, the main difficulty lies 
in computing the gravitational field of a given mass distribution accurately and 
efficiently. Any gravodynamical code must include at least this much, including 
those used for studying planet and cosmological structure formation. 

Many algorithms have been devised for calculating gravitational fiel ds of iso- 
lated galaxies. The simplest is direct particle-particle interactions (e.g. lAarsethl 
1963), but this scales as 0(N 2 ), N being the number of particles, making it 
impractical for other than special applica tions. Particle- Mesh (PM/grid) codes 
(|Miller fc Prendergastlll96^ : lMillerlll976t ) solve the potential on a grid by bin- 
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ning the mass distribution. Several geometries are possible, including carte- 
sian, cylindrical and spherical. Grid codes are very efficient, scaling largely 
with the number of grid cells, N g , as ~ 0(N g ) (Sel lwoodlll997r ), Their main 
drawbacks are their inflexible geometry and their trade-off between volume 
and spatial resolution, but hy brid codes can alleviate many of these problems 
fe.g.lFudl!999l:ISellwoodll2nn.'ir). Adaptive Me sh Refinement (AMR) codes (e.g. 
iBrvan k Normanlll995l : iKravtsov et al.lll997f ). on the other hand, let the grid 
evolve dynamically t o increase the re solution in dense regions. Tree codes (e.g. 
iBarnes k Hudll986l : lHernouisdll987t ) instead group particles by location, com- 
puting direct forces for nearby particles and a few multipoles for the more dis- 
tant ones. These codes can be vectorized effic iently and scale as Q (N dog N) 
(|DubinskJl99ilStadel200ll : ISpringel et aLlfcOOlT ) while lDehnenl (l200ft > presented 
an 0{N) extension. Tree codes have been widely implem ented and several are 
publicly available. Self-Consistent Field (SCF) cod es (e.g. IClutton-Brocld 119721 : 
iHernquist k Ostrikerlll992l : IE arn k Sellwoodlll995l l expand the density in a set 
of orthogonal basis functions, from which the potential can then be computed. 
SCF codes scale as 0{N) and do not require softening because truncating the 
expansion is enough to suppress small scale noise. However they require that 
the basis set be chosen with care or that a large number of basis f unctions be 
included (which re- introduces small scale noise). IWeinberd (Il99fih solves this 
problem by making the basis set adaptive. 



2.1. Code testing 

The first requirement of any TV-body code used to study galaxy evolution is that 
it be collisionless. R elaxation rates of conserved quantities are useful diagnostics 
for th is purpose (e.g.lHohllll973l:lHernouist k Barnesll9"90l : lHernauist k Ostrikerl 
119921 : IValenzuela k KlvpinlfcOO.t t . 

Testing the gravodynamical part of a code requires more effort. Compari- 
son with the limited number of exact analytic results know n provide strin gent 
code tests; s uch systems (in 2-D) include the Kalnajs disk dKalnaislll972T ). the 
Mestel dis k (IZand Il976h, the i sochrone disk ( Kalnaisl Il978t ) and the power- 



law disks ((Evans k R ead 1998a. bj). Several co de tests using these predictions 
have been carried out. lEarn k Sellwoodl (|l995l ) compared the eigenfrequencies 
of instabilities in isochrone disks with simulations using an SCF and a PM 
method. They were able to reproduce the predicted values to within 5% with 
just 15K particles using the SCF method whereas softening c aused the PM code 
to nev er quite converge to the predicted value. Meanwhile, ISellwood k Evansl 
( 2001) presented TV-body examples of the power-law disks including reproduc- 



ing a challenging case of a perfectly stable disk. They reported that aliasing 
rendered SCF codes ill-suited to these disks. Insta bilities of s pheri cal systems 
have also been used for code testing. For example Weinber 3 rtl99fiT) tested his 
adaptive S CF code on the inst abilities of spherical generalized polytropes inves- 
tigated bv IBarnes et al.l ( 19861 ) and found generally good agreement. Analytic 



solutions of no n-equilibrium evo lution, which include the solution of 1-D plane 
wave collapse dZel'Dovichl ll970T) and of spherical in fall in an expanding uni- 
verse dFillmore k OoldreicUli98l lBertschingerlll985r) have also be used to test 



N-bodv codes (|Kravtsov et aljll997l : iDave et al.lll997fi . 
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Different codes can be tested also by direct comparison. Ilnaeaki et al.l 

( 1984D compared bar formation between an iV-b ody model and a direct numeri- 
cal in tegration of the CBE and Poisson equation (|Nishida et al.ll98ll : IWatanabe et al.l 
Il98lh . The two simulations matched each other to better than 2% in bar ampli- 
tude well into the non-linear regime and eventual discrepancies were due to the 
inability of the CBE code to handle large gradients in the distribution function 
(DF). ' 



3. Cusp evolution 



A prediction of cold dark matter (CDM) cosmology is that dark m atter halos are 
cusped, with densities p y at small radii and 1 < < 1.5 ( Navarro et alJ 
119971 : iMoore et all Il998l : I.Ting fc Sutol 120001 : iPower et al.l l2003T t . Several argu- 
ments have been advanced against the presence of such cusps in real galax- 
ies, incl uding detailed rotation curve fits of dwarf and low surface brightness 
galax ies (|Blais-Ouellette et al.ll200l|:lde Blok et alJl200ll:lMatthews fc Gallagherl 
2002) , bar gas flows (IWeiner et al.ll200ll : iPerez et al.ll2004l ) and pattern speeds 
(|Debattista Sellwoodl ll998. 2000;, see fllH b elow) and the microlensing optical 
depth and gas dynamics in the Milky Way (|Binnev et alJ 120001 : iBissantz et alJ 
120031 ). 

Sev eral ways to erase cusps have been considered including new dark matter 
phys i cs dSpergel fc Steinhardtl2000l : lPeeblesl2000l:lGoodmarJ2000ljKaDlinghat et al l 
200C : ICenl200lh . feedback from star formation ([Navarro et alll996l : lGnedin fe Zhad 



2002) . or an initial power spectrum with d ec reased power on sma ll scales ((Hogan Sz Dalcantonl 



2000: ICohn et al.ll2000l : bode et al.ll200ll ). iBinnev et all (l200lh suggested that 



bars in young galaxies were able to torque up cusps and expel them, while 
iHernauist k, Weinberg! (|l992f ) had found that such torques can reduce the den- 
sity of a spheroid by a factor of ~ 100 out to ~ 0.3ob (where a# is the 
bar's semi-major axis) . IWeinberg fc Katd (l200l hereafter WK02), following 
IHernauist fc Weinberg) (|l992l ). presented perturbation theory calculations and 
simulations of imposed rigidly-rotating non-slowing (IRRNS) bars to argue that 
cusps are erased in a few bar rotations (~ 10 8 years). Their bars were large (~ 10 
kpc) , but they argued for a scenario in which large primordial bars destroy cusps 
while the current generation of smaller bars formed later. Large reductions of 
dark halo densities were not seen in self-consistent c uspy halo simulations (e.g. 
lO'Neill fc Dubinskil2003l : IValenzuela fc KlvDinl2003l ) : WK02 suggested that this 
was because the resonant dynamics responsible for cusp removal are very sen- 
sitive to numerical noise. They found that they needed > 10 6 for their SCF 
code and argued that even larger TV" would be needed for grid, tree or direct 
codes. 

The scenario of WK02 makes three main claims: (1) that bars cause a 
decrease in cusp density and can destroy cusps if sufficiently large and strong 
(2) that such bars formed via interactions at high redshift and (3) that N needs 
to be large in order that the relevant phase space is adequately covered and that 
orbit diffusion does not destroy the resonant dynamics causing cusp destruction. 
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3.1. When do bars destroy cusps? 

ISellwoocO (|2003l . hereafter S03) presented simulations of the same IRRNS bar as 
was used by WK02 and found cusp destruction occurred in a runaway process. 
The inclusion of the I = 1 spherical harmonic terms in the potential had a large 
effect on the evolution; in their absence, cusp removal took ^5 — 6 times longer 
than w hen they were present. For the same IRRNS bar, fMcMill anfc Dehnenl 
( 2005, hereafter MD05) found that the cusp moves off-center by as much as 30% 
of as. This centering instability gives the appearance that the cusp has been 
erased when halo density is measured by spherically averaging about the origin 
which, they argued, caused WK02 to over-estimate a bar's ability to erase a cusp. 
Evidence for lopsidedness in the simulation of WK02 can be inferred from its 
asymmetric bar-induced halo wake (their figure 2). When MD05 suppressed this 
purely numerical instability they still found that the cusp is destroyed, although 
after a much longer time. The possibi lity that an offset cusp is being co nfused for 
cusp destruction was investigated by iHollev-Bockelmann et al. (2005, hereafter 
HB05). Their self-consistent simulations formed bars in a disk of particles, thus 
damping any centering instability. While confirming that some of the evolution 
seen by WK02 was due to the centering instability, they found that when odd I 
terms were excluded in their simulations that the damage to the cusp was not 
significantly diminished. 

Thus these idealized bars at least are able to destroy cusps. What about 
more realistic bars? S03 showed that when, instead of assuming a fixed pattern 
speed, f2 p , he allowed it to decrease such that total angular momentum is con- 
served assuming a constant moment of inertia (an IRRS bar), that the cusp was 
damaged significantly less. This happens because an IRRNS bar must transfer 
more angular momentum to the halo than it can plausibly have in order to de- 
stroy the cusp. Shorter, more realistically sized bars were also ineffective bar 
destroyers. Similarly MD05 found that an IRRS bar is unable to destroy a cusp 
in a Hubble time. S03 also presented self-consistent simulations in which a disk 
of particles was grown inside a halo. In these cases, besides forming smaller bars 
which slow, cusp erasure was also inhibited by the increase in halo density as the 
disk grew and again as the central density of the disk increased because of angu- 
lar momentum transport outwards. HB05 instead found that the large relative 
angular momentum gained by the inner halo in their self-consistent simulations 
flattened the cusp, but only out to 900 pc scaled to the Milky Way. 



3.2. Do bars get as large as needed? 

Thus all self-consistent simulations find that bars of sizes typical of those ob- 
served are unable to remove cusps on scales of several kpc. Do b ars ever get to 
be large enough to do so? Bars that form via disk instabilities (|Toomrelil98lh 



generally extend to roughly the radius of the rotation curve turnover. However, 
externally triggered bars may be substantially larger. HB05 showed an example 
of such a bar; scaled to the Milky Way, external triggering produced a bar of 
as = 12 kpc, as opposed to ~ 4.5 kpc via the bar instability. Whether the re- 
quired triggering actually occurs is a question for hierarchical simulations while 
direct observations at high redshift should establish whether bars get to be as 
large as 10 kpc. 
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3.3. Is the evolution in simulations compromised by too small iV? 

I Weinberg Kat3 ((2005h provided several reasons for the need of large N . They 
focused especially on the inner Lindblad resonance (ILR) which always extends 
down to the cusp. One possible problem they identified is two-body scattering 
which leads to particles executing a random walk and therefore lingering in 
resonance for significantly les s than they would otherwise. In the terminology of 
iTremaine &: Weinberd (|1984). scattering causes particles to traverse the ILR in 
the fast rather than the slow regime. Since resonant torques scale as mp in the 
fast regime and as yjmp in the slow regime (where mp is the fractional mass 
of the perturber) two-body scattering causes the friction to be substantially 
reduced. For a typical bar, they argued that N > 10 8 within the virial radius 
is required. A second limiting factor they identified is phase space coverage. 
Whether a resonant particle gains or loses angular momentum depends on its 
phase. If the phase space of the resonance is not sufficiently sampled by particles 
then the correct ensemble average is not attained. 

HB05 found no difference between their tidally perturbed simulations with 
N = 5.5M and N = 11M, while N = 1.1M resulted in a weaker bar and less 
friction than in the other two cases. This, they argued, proved that the N = 
1.1M simulation su ffered from too much orbital diffusion to follow resonances 
correctly. However, ISellwoodl (|2005h argued plausibly that this N dependence 
was due to the bar being weaker, and not because friction depends on N. The 
weaker bar in the N = 1.1M simulation was caused by the unavoidably larger 
m = 2 seed amplitude in the initial disk which, once swing-amplified, interfered 
destructively with the distortion induced by the externally applied tidal field. 

The simulations of S03 (using both a PM and an SCF code) and MD05 
( using a tree code) are instructive because they studied the same physical system: 
a iHernauistl (|l990l ) halo with a IRRNS bar of cle = 0.7r s (r s being the halo 
scale radius). Both studies found that the onset of the runaway cusp destruction 
depended on N and occurs earlier with decreasing N. S03 found good agreement 
between evolution using the PM code and the SCF code; additionally he showed 
that using only I = and 2 terms was sufficient for the evolution, with little 
change when larger even values of I were included. On the other hand, MD05 
forced symmetry about the origin to suppress the centering instability. Thus the 
two studies were very similar; but while their results are in qualitative agreement, 
there are surprising quantitative differences. S03 found that the onset of the 
runaway appears to be converging by N = 3M at ~ 130 bar rotations (his 
figure 2b) but for the same N MD05 saw no evidence for convergence with 
runaway at ~ 210 bar rotations (their figure 5b). Yet both agreed that the 
radius containing 1% of the total mass is driven from ~ 0.1r s to ~ 0.5r s . Unless 
the difference is due to the initial conditions in some worrisome way, noise would 
seem to responsible for these differences. A more careful comparison between 
thes e different simulations seems p articularly worthwhile. 

ISellwood Sz Debattistal (|2005h report an entirely different test of whether 
scattering overwhelms bar evolution. They presented simulations (using a PM 
code) in which they perturbed the system into a metastable state in which 
J7 p was nearly constant (more on this in iJ4.|) . Their metastable state resulted 
because £l p was driven up such that the principal resonances are trapped in 
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shallow local minima of the DF. Systems persisted in this metastable state for 
~ 5 Gyr, which would not have been possible if orbit scattering had been strong. 

Unlike two-body scattering, phase space coverage can be investigated by 
means of test particle simulations, whe re the sel f- gravit ating halo response is 
not included. Using such experiments, ISellwoodl ((20051 1 found that N ~ 1M 
was sufficient for the fJ p evolution of a low mass bar Mb ar = 0.005Mh a i o to 
converge, and a factor of ~ 100 less particles were needed for Mb ar = 0.02Mh a i o - 
When self-gravity was introduced (thus adding scattering), 1M particles were 
needed for the Mb ar = 0.02Afh a i o case. As in S03, he argued that WK02 needed 
large N because of the extremely difficult nature of their experiments with a 
fixed-amplitude bar rotating at a fixed O p . When O p is a function of time, the 
resonances are broadened to orbits near resonance, which makes phase space 
coverage a less stringent constraint on N. Weak bars in particular may be prone 
to such difficulties, but it is strong bars which are most interesting. 

In summary, there may be some evidence that noise is compromising some 
if not all of the evolution in some cases, in ways not yet fully understood. But 
it is clear that only quite large bars can erase cusps to scales larger than 1 kpc. 



4. Evolution of pattern speeds 

The pattern speed of a bar is usually parametrized by the ratio 7Z = D^/as, 
where Dl is the corotation radius, at which the gr avitational and cen trifugal 
forces cancel. A self-consistent bar must have K > 1 (IContoDoulosl H980 N l. A bar 
is termed fast when 1.0 <1Z< 1.4 and slow otherwise. This definition does not 
distinguish fast from slow in terms of Q p alone: a bar in a galaxy with rotation 
velocity of 200 km s~ x is slow at O p = 100 km s~ x kpc^ 1 if cib = 1 kpc, but 
fast at flp = 20 km s _1 kpc -1 if clb = 10 kpc. O bservational evidence points 
to fast bars, both i n early-type barred g a laxies (iMerrifield fc Kuiikenl Il995l: 
Gerssen et all Il999l: iDebattista etafl 120021 : lAguerri et alJ 120031: iGerssen et all 



200.1: iDebattista fc Williamsl2004h an d in late-types (iLindblad et al.ll99fil : lLindblad fc Kristenl 
199ft IWeiner et al.ll200ll: iPerez et al.ll2004h . 

rTre^mdn^^^Weinljera ((1984 ) d"eveloped the perturbation theory of dynami- 
cal friction for perturbers in spheroidal systems, showing that friction arises near 
resonances, when mO p = kO, r + lQ^, where O r and are the radial and angular 
frequencies, resp ectively fsee IWeinberdl2005l . for a time-dependent treatment). 
IWeinbergl ( 1985T ) applied this theory to a bar rotating in a dark halo and found 



that the bar is braked such that Dl 3> clb unless (1) angular momentum is 
added to the bar, (2) the bar is weak or (3) the halo has low mass. 

The transfer of angular momentum fro m disk to spher o id (bulge or halo) 
was reported in several early sim ulations ((Sellwoodl Il980l : iLittle h Carlbera 
ll99ll : lHernauist k, Weinberdll992f ). Fully self-consistent, 3-D simulations of bar- 
unstable dis ks embedded in dark halos were presented bv lDebattista k, Sellwoodl 
((19981 . l2000l . hereafter DS00) who found that disks needed to be the dominant 



has also been found in other simulations since then ( 


O'Neill fc Dubinskil 12003 ; 


Hollev-Bockelmann et al.ll2005l:lBerentzen et al.ll2005l: 


Martinez- Valpuesta et al.| 


2005). lO'Neill k Dubinskil ((2003) included a comparison with DS00 and found 



generally good agreement. 
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T he constraint of DSOO severe ly limits CDM cusps to be present in dark 
halos. IValenzuela Sz Klvpinl (|2003l . hereafter VK03) presented simulations of 
a CDM Milky Way, with halo concentration c ~ 15. These simulations were 
evolved on an AMR code with a maximum refinement corresponding to a res- 
olution of 20-40 pc. They found that the bar which formed remained fast for 
~ 4 — 5 Gyr, thereafter slowing to 1Z = 1.7. They concluded that the slow bars in 
the simulations of DSOO were an artifact of low resolution. They also argued that 
their bars remained modest in length (1-2 scale- lengths) but became too long in 
lower resolution simulations, which they concluded again pr oved the limitations 
of low resolution simulations. ISellwood fc Debattistal (|2005l repeated one of the 
simulation s of VK03 usin g the same initial conditions but evolved on a hybrid 
grid code ([Sellwoodl 12003) with a fixed softening comparable to that of VK03. 
They found instead that the bar reached 1Z = 2 within 4 Gyr. They accounted 
for the result of VK03 by noting that as the bar formed, the central density 
increased (see fl5~j) . Thus in an AMR code there is a tendency for the spatial 
resolution to increase; the disk in the model of VK03 being rather thin, this 
led to enhanced forces and a corresponding artificial increase in Q p . Once this 
happened, the bar found itself with its principal resonances at local phase space 
minima previously generated by the forming bar. Although angular momentum 
is exchanged at resonances, the sign of the torque depends on the phase space 
gradient of the DF at the resonances; in the absence of a gradient no friction 
is possible. Thus the numerics induced a metastable state. This state proved 
fragile to realistic perturbations and is not likely to last long in nature, but in 
the quieter environment of an isolated iV-body system it can persist for many 
Gyrs, enough to fully account for the behavior found by VK03. 

In a series of papers, Athanassoula present ed several arguments against 
the conclusion of DSOO. In lAthanassoulal (2003) she compared the evolution 
of two halo-dominated systems, MQ2 and MHH2, with nearly identical disk 
and halo rotation curves. However, their velocity dispersions were different, 
being larger in MHH2 because its halo extended to larger radii. As a result O p 
decreased significantly in MQ2 but hardly at all in MHH2. She concluded that 
this "argues against a link between relative halo content and bar slowdown..." 
This claim however is contradicted by her own simulations. The confusion arises 
from her relying on the change in J7 p to constrain the halo. Not only is this 
unobservable, but if we compute 1Z for her models, we find 1.4 < 1Z < 1.7 for 
run MQ2 and 1Z ~ 3 for MHH2, i.e. both these halo-dominated systems are 
slower than observed. Far from being in disagreement with the results of DSOO, 
her simulations support them. She also showed that weak bars are less able to 
drive angular momentum exchanges and may remain fast; since observational 
mea surements have only been obtaine d for strong bars, this is not worrisome. 

lAthanassoula &: M isiriotis (2002) argued that bar lengths are difficult to 
measure in simulations making a comparison with corotation difficult. Bar 
lengths are certainly not always easy to measure, but it is still possible to define 
values which straddle the real value; for example DSOO had cases in which the 
uncertainty in a n was as large as A or/or ~ 0.3, not substantially different than 
in observations ( Debattistal 12003) . But these are not Gaussian errors and the 
probability of as falling outside the given range is practically zero. Moreover, 
the same problem afflicts observations; therefore measurement uncertainties in 
olb are a nuisance but not a repudiation of the constraint. 
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lAthanassoulal (12002ft showed that loss of angular momentum from the bar 
leads to a growing bar, as first suggested by DSOO. VK03 pointed out that bars 
become excessively long in the presence of strong friction. Even though angular 
momentum redistribut ion leads to larger disk scale-lengths, bars extending ^ 10 
kpc are not common (|Erwinl l2005). Strong bar growth does not seem to have 
occurred through the history of the current generation of bars and presumably, 
neither has strong friction. 

Thus to date no well-resolved simulation has p rovided a valid counter- 
example to the claim bv iDebattista fe Sellwoodl £l998) that dense halos cannot 
support fast strong bars. If anything, recent simulations have lent support to it. 



The secular evolution of disk densities 



The excellent recent review bv lKormendv Sz Kennicuttl ( 2004f ) presents the ob 



servational evidence for pseudo-bulge formation via secular evolution and dis- 
cusses some of the older A-body results in that field. Here I review recent 
developments not covered by those authors. 

5.1. Disk Profile Evolution 

Ba r formation is accompanied by a rearrangement of disk material as first shown 
bv iHohll ( 197ll ). Generally a nearly double-exponential profile develops, with a 



smaller central scale- length than the initial and a larger one further out. In this 
respect these profiles resemble bulge+disk profiles and comparisons with obser- 
vatio ns show that these profiles are reasonable appro x imations to obser v ed pro - 
files (|Debattista et al.ll2004l : lAvila-Reese et al.ll2005l ). IDebattista et alJ (|2005bh 



showed that the degree by which the profile changes depends on the initial disk 
temperature Q. In hot disks {Q ~ 2) little angular momentum needs to be shed 
and the azimuthally-averaged density profile is practically unchanged. Thus the 
distribution of disk scale-lengths depends not only on the initial angular mo- 
mentum of the baryons (and presumably that of the dark halo) but also the 
disk temperature. Angular momentum redistribution continues also after the 
bar forms, especially to the halo. This leads to a further increase in the cen- 
tral density of the disk even when the evolution is collisionless. This may be 
sufficient to render an initi ally halo-dominated system into one dominated by 
baryo ns in the inner parts ( Debattista h Sellwoodl l2000l : IValenzuela &: Klvpinl 



12003ft . 

The angular momentum lost by the bar as it forms may be transported out 
to larg e radii via a resonant co upling between the bar an d spir als (IDebattista et al.l 
l2005bli of the kind found bv iMasset Taggerl (|l997fi and iRautiainen Said 



1999). IDebattista et al.l (l2005bh show that this transport leads to breaks in 



the density distribution which, when viewed edge-on, are indistinguishable from 
those observed in real galaxies (lPohlenll2002l N l. 



5.2. Vertical Evolution: Peanut-Shaped Bulges 

A-body simulations of the vertical evolution of disk galaxies have concentrated 
on box- or pean ut- (B/P-) shaped b ulges which are present in some 45% of 
edge-on galaxies ( Liitticke et al.l l2000). Simulations have shown that these form 
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via the secular evolution of bars ( Combes k, Sanderslll98ll ) , eith er through reso- 
nant scattering or t h rough bending (aka buckl i ng) instabilities (lPfennigerlll984l : 
ICombes et all Il990l : iPfenniger Friedlil Il99ll : IRaha et alJ Il99lh . Although it 
is often thought that a peanut requires that a bulge was built by secular pro- 
cesses, simulations show that peanuts can also form when the initial conditions 
include a bulge, as would happen if bulges form through m ergers at high redshift 
( Athanassoula fc Misiriotisll2002l : iDebattista et aLlfcOOBaT l. 

Obs ervations seeking to establish the connection between B/P-shaped bulges 
and bars dKuiiken Sz Merrifieldll995l : iMerrifield &: Kuiikenll999l : iBureau &: Freeman! 
1999; IChung fc Bureaul 12004) . by looking for evidence of bars in edge-on B/P- 

bulged systems, have benefited f rom comparisons with the edge-on s te llar veloc- 

ity dis tributions of iV-body bars ( Bureau &i Athanassoulall999l . l200'5l ). IBureau &: Athanassoulal 
(2005) characterized the signature of an edge-on bar as having (1) a Freeman 
type II profile, (2) a rotation curve with a local maximum interior to its flat 
part (3) a broad velocity dispersion profile with a plateau at intermediate radii 
(4) a correlation between velocity and the third -order Gauss- Hermite moment 
h 3 ((Gerhard! I"l993: van der Marel Sz Franxlll9 93). A diagnostic of B/P-shaped 
bulges in face-on galaxi es was developed, and tes ted on high mass and force res- 
olution simulations, bv IDebattista et "all ( 2005al ). Vertical velocity dispersions 
constitute a poor diagnostic because they depend on the local surface density. 
Instead, their diagnostic is based on the fact that peanut shapes are associated 
with a flat density distribution in the vertical direction. The kinematic sig- 
nature corresponding to such a distribution is a minimum in the fourth-order 

Gauss-Hermite moment /14. 

Th e buckling instability itself has also been studied with simulations. IDebattista et all 
(2005b) showed that an otherwise vertically-stable bar is destabilized when it 
slow s and grows. Moreover, the insta bility can occur more than once fo r a given 
bar ((Martinez- Valouesta et all 120051 ) . Finally, after IRaha et alJ (|l99lh showed 
that bucklin g weakens bars, it has of ten been assumed that bars are destroyed 
■ IDebattista et aD ( 2005bh presented a series of high force and mass 



by buckling 

resolution simulations demonstrating that this is not the case. 



5.3. Spirals 



There is broad agreement that spirals constitute density waves (|Lin fc ShiJl964h 



At least three different dynamic al mechanism s have been propos ed for excit- 
ing th em: swing-amplification ( Toomrelll98lh . long-lived modes (|Bertin et al.l 



1989aj0) and rec urrent instabilities seeded by features in the angular momen - 
tum distribution (jSellwood fc Linlll989l : ISellwood fc K Ih^ll99ll : ISerlwoo 



All three are still v iable and not much new A -body results have been obtained 
in recent years, but lSellwood fc Binnevl |2002) used simulations to demonstrate 



that spirals cause a considerable radial shuffling of mass at nearly fixed angular 
momentum distribution. This happens at corotation and is not accompanied by 
substantial heating — stars on nearly circular orbits can be scattered onto other 
nearly circular orbits. For example, they estimate that a star born at the solar 
radius can be scattered nearly uniformly within AR = ±4 kpc. Thus the idea 
of a Galactic habitable zone becomes somewhat suspect, as does the need for 
infall to maintain the metallicity distribution observed at the solar circle. 
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6. Bar destruction by CMCs 

Central massive concentrations (CMCs), whether supermassive black holes (hard 
CMCs) or gas condensations several hundred parsecs in size (soft CMCs), could 
destroy bars. Early studies of this phenomenon were ins pired by the simi- 
lar w ork for slowly-rotating triaxial elliptical galaxies (e.g. iGerhard Binnevl 
1985), where the loss of triaxiality results from the destruction of box orbits by 
scattering off the CMC. Such scattering may not be efficient in bars, since the 
main bar-s upporting orb i ts are loops which avoid the center. IHasan Norman] 
(1990) and IHasan et al.1 (|l993T l argued that, when a CMC in a barred galaxy 
grew sufficiently massive, it quickly destroys bar-supporting orbits by driving 
an ILR, around which orbits are unstable, to large radii. The more centrally 
concentrated the CMC grew, the further out was the ILR and therefore the more 
disruptive it was. How massive the CMC needed to be required iV-body simula- 
tions to establish. iFriedli fc Ben^ ((l993h modeled gas and stars with PM+SPH 



simulations and found that gas i nflows destroyed bars when 2% of the baryonic 
mass ended up in a hard CMC. iNorman et alJ (1996) pursued collisionless 2-D 



and 3-D simulations, with CMCs grown by slowly contracting a massive compo- 
nent. They needed a CMC of mass 5% of disk mass, Md, to destroy bars. Then 
bar destruction was rapid an d led to a bulge-like spheroid. 

IShen fc SellwoodI (120041 . hereafter SS04) presented a series of high-quality 
simulations including high mass, force and time resolution and found that bars 
are more robust than previously thought. They varied the growth rates, com- 
pactness and mass of the CMCs and considered both weak and strong bars 
within which they grew CMCs at fixed compactness. A rigid halo with a large 
core radius was also included. They obtained a fast decay while a CMC was 
growing followed by a more gradual decay once the CMC reached its full mass. 
The time over which the bar was grown proved unimportant, with only the final 
mass and compactness mattering. Hard CMCs cause more damage than soft 
ones, needing 4 — 5%Md and > 10%Md, respectively, to destroy bars. Their 
tests showed that time steps need to be as small as 10 -4 of a dynamical time 
in order that more rapid but incorrect bar destruction is avoided. They inter- 
preted the two-phase bar destruction as scattering of low energy bar-supporting 
x\ orbits during the CMC growth phase and continued gradual global structural 
adjustment, which further destroyed high energy x\ o rbits, thereafter. They 
predicted that massive ha los, which lead to bar growth ( Debattista SellwoodI 



2000: lAlhanassoulal l2003f ) render bars even more difficult to destroy; this pre- 
diction, as well as the two-phase bar w eakening, was confirmed in live-halo sim- 
ulat ions bv lAtha nassoul a et alJ ((20051) . 

iBournaud Sz Combed ([2002 ) and lBournaud et alJ ((20051 ) have argued for a 
radically different picture. They simulated gas accretion and noted that bars 
were destroyed and reformed 3 or 4 times over a Hubble time. The amount of 
gas required to destroy bars is not, however, wholly consistent in these simu- 
lat ions: a system with ~ 7% total gas mass fraction lost its bar within 2 Gyr 
in IBournaud et al. (2005) whereas a s ystem with three times more gas main- 
tained its bar in IBournaud &: Combed (|2002f ). Possibly star-formation, included 
in the later simulations, someh ow quenched the infal l onto the center. SS04 
hinted that the t i meste p used bv IBournaud h Combed ( 20021 ) was too large but 



IBournaud et alJ ( 2005h reported that using a timestep 0.125x their standard 
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value (and close to that adv ocated by SS04) still led to recurrent bar destruc- 
tion. iBournaud et al.l (|2005h argued that their results are correct and that other 
studies had erred in mimicking gas accretion by simply growing a massive ob- 
ject because this neglected the important effects of angular momentum trans- 
port from gas to the b ar. This is, perhaps, consistent with the simulations of 
iBerentzen et al.l (|l998f ) (who were able to destroy bars b y gas inflow lead i ng to a 
CMC of mass fraction just 1.6%) and the earl ier ones oflFriedli fc Bend (Jl993). 
On the other hand, the fully live simulations of lDebattista et alJ ((2005bh only de- 
stroyed bars when soft CMCs reached ~ 20%M<j, in good agree ment with SS04. 
One p ossibly important difference is that the simulations of IBournaud et al.l 
(2005) included a rigid halo, which prevents it fr om accepting angular momen- 
tum frOTi_the_j3a£ J _w^^ (2005b) had live halos. 

iHozumi fc Hernquistl (l2005T l. using a 2-D SCF code, also concluded that 
hard CMCs ca n destroy bars with smaller masses, 0.5%Md. SS04 (see also 
ISellwoodl 120021 ) speculated that low order SCF expansions may not be able 
to simultaneously maintain a system axisymmetric near the center and non- 
axisymmetric further out. 



7. Multiple patterns 

An emerging field in the past few years has been galaxies with multiple pat- 
terns. These are challenging to study because traditional tools such as surfaces- 
of-section are no longer viable. TV-body simulations, therefore, are vital for 
studying systems such as bars in triaxial halos and bars within bars. 



7.1. Bars in triaxial halos 

CD M predicts that dark matter halos are triaxial Jb arnes _fc EfstathioT] Il987l : 
iFrenk et al.lll988l : iDubinski fc Carlbergll99ll : I.Ting fc Sutdl2002T ). The omdensa- 
tion of baryons inside triaxial halos drives the m to rounder shapes, but systems 
do no t generally become wholly ax isvmmetric (lDubinskilll994l:lK'azantzidis et alJ 
Using rigidly-rotating bars. lEl-Zant &: Shlosmanl (|2002r ) computed the Li- 
apunov exponents of orbits and showed that chaos quickly dom inates the evolu- 
tion w hen halos are triaxial and cuspy. iV-body simulations bv IBerentzen et"aTI 
(2005) indeed show that bars are destroyed in cases where the triaxiality in the 
potential is as small as c/o ~ 0.9 and the halo is cuspy. One way in which this 
fate can be avoided is for the bar to alter the shape of the inner halo. At present 
it is not clear which systems can accomplish this and which cannot. A better 
understanding of when bars are destroyed in triaxial halos could possibly lead 
to an important new constraint on the shapes and profiles of dark matter halos. 



7.2. Bars within bars 

While observations of largely gas-free early-type galaxi e s find an abundance 
of nuclear bars within large scale bars ( Erwin Sparkel (|2002l ) found them in 
~ 30% of barred SO-Sa galaxies), simulating them has proved difficult. Moreover, 
it is only recently that direct observational evidence for kin ematically decouple d 
primary and nuclear bars in one system has been obtained ( Corsini et al.l l2003). 
Therefore their dynamics have been poorly understood, in spite of the fact that 
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they have been postulated to be an import ant mechanism for d riving gas to 

small radii to feed supermassive black holes ( Shlosman et al.lll989t ). 

Most numerical studies have required gas to form second ary bars ( Friedli Sz Martinet! 
119931 : IShlosman k Helledl2002l : lEnglmaier fc Shlosmanll2004h , but their presence 
in gas-poor early-type galaxies suggests that gas is not the main ingredient for 
formi ng secondary bars. Stellar counter-rot a tion can lead to counte r-rotating 
bars (|Sellwood fc Merrittlll99l : lFriedlilll99fiT l. iRautiainen et al.l {2002) were the 
first to succeed in producing collisionless iV-body simulations with both bars 
rotating in the same sense. Their secondary bars, which were vaguely spiral-like 
and possibly hollow, rotated faster than the primary bars and survived for sev- 
eral Gyrs. Debattista & Shen (2006, in preparation, see also Shen & Debattista 
in these proceedings) present further examples and explore the mutual evolution 
of the two bars. Now that iV-body simulations can achieve the high force and 
mass resolution needed to form self-consistent nested bar systems it is hoped 
that progress in understan ding these systems wil l be more rapid than in the 30 
years since their discovery ( de Vaucouleursiri975l b 



8. Desiderata for the future 

The collisionless simulation of isolated galaxies is an endeavor over thirty years 
old. The subject is still rich, with several open problems, and continues to be 
very active. Algorithmically, if not conceptually, it is the simplest problem that 
can be considered. Various gravity solvers for its study are available (which are 
used also in other areas of astronomy). Despite much progress, the degree of 
disagreement in the field, as described above, is surprising. The way to progress 
from this point is to compare directly the results of different codes with eac h 
other, as has been done in other fields (e.g. iKang et al. l ll994l : lFrenk" et al. 1999). 
Unfortunately the number of such tests in galaxy evolution has been small. 
Therefore a detailed comparison between many different codes would be very 
valuable at this time. Ideally this would involve several implementations of the 
same code type to establish behaviors in the different types. The actual tests 
to be performed should include systems in which an analytic result is known 
(useful to establish values of numerical parameters which are optimal) as well 
as systems for which the result is not known in advance. Furth e rmore , the N- 
body tests of the type recently proposed by IWeinberg fc Katzl J2005) provide 
a challenging and useful basis for assessing the effect of noise on simulations. 
Such a comparison will give the community greater confidence that we are able 
to model correctly the most basic level of galaxy evolution. 
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